Controle de processos βAR(1)

# Temas para os gráficos

theme.base <- theme_minimal(base_size = 11) +
  theme(
    axis.text = element_text(size = 8),
    plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 9),
    plot.caption = element_text(size = 8),
    axis.title = element_text(size = 8),
    legend.title = element_text(size = 8),
    panel.grid.major = element_line(colour = "grey90", linewidth = 0.5),
    panel.grid.minor = element_line(colour = adjustcolor("grey90", alpha.f = 0.4), linewidth = 0.5),
    panel.border = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.line.x = element_line(colour = "grey"),
    axis.line.y = element_line(colour = "grey"),
  )

theme.no_legend <- theme(legend.position = "none")

theme.no_grid <-  theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank()
)

theme.no_axis <- theme(
  axis.line.x = element_blank(),
  axis.line.y = element_blank()
)

theme.no_axis_title <- theme(
  axis.title.x = element_blank(),
  axis.title.y = element_blank()
)

# Theme for timeseries
theme.ts <- theme.base +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  ) +
  theme.no_legend

plotly.base <- function(p) {
  p %>%
    layout(
      margin = list(b = 60, t = 80)
    ) %>%
    config(mathjax = 'cdn')
}

# Função para salvar e carregar resultados
cache_dados <- function(chave, funcao_geradora) {
  c <- paste(chave, "rds", sep = ".")
  cache <- file.exists(c)
  
  if (!cache) {
    resultado <- funcao_geradora()
    saveRDS(resultado, c)
  } else {
    resultado <- readRDS(c)
  }
  
  resultado
}

Introdução

βAR(1) é um modelo de séries temporais autorregressiva de ordem 1 com distribuição Beta – onde a variável de interesse está restrita no intervalo (0, 1) –, este modelo foi proposto por Rocha e Cribari-Neto:

• Beta autoregressive moving average models. 2009. DOI: https://doi.org/10.1007/s11749-008-0112-z
• Erratum to: Beta autoregressive moving average models. 2017. DOI: https://doi.org/10.1007/s11749-017-0528-4

O objetivo do presente trabalho é aplicar métodos de controle de processos βAR(1) para detectar mudanças no processo. Para isso, utilizamos o fato de que, quando modelada corretamente, os resíduos de uma série temporal, como a βAR(1), são independentes e normalmente distribuídos com média zero e variância constante.

Desta forma, quando o processo sofre uma mudança, espera-se que os resíduos não sejam mais independentes e que a média e a variância dos resíduos sejam diferentes. Assim, podemos utilizar métodos de controle de processos para detectar essas mudanças.

Modelo

Olhando sob a perspectiva de teste de hipóteses, temos:

• \(H_0\): amostra inicial de tamanho 100, com \(\Phi_0 = 0.2\). Indicando a Fase I da aplicação do controle de processos.
• \(H_1\): amostra subsequente de tamanho 200, com \(\Phi_1 = 0.2, 0.3, \ldots, 0.6\). Indicando a Fase II da aplicação do controle de processos.

Além disso, utilizaremos um nível de significância de 0.05 para os testes de hipóteses, o que nos dá um quantil de 1.96 para a distribuição normal padrão.

A simulação dos dados será feita a partir da biblioteca BTSR: Bounded Time Series Regression.

nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2
alpha_teste <- 0.05
quantil_teste <- qnorm(1 - alpha_teste / 2)

coeficientes <- function(phi) {
  # βARMA: the model from Rocha and Cribari-Neto (2009, 2017)
  #        is obtained by setting `coefs$d = 0`
  #        and `d = FALSE` and `error.scale = 1` (predictive scale)
  list(alpha = 0, nu = 20, phi = phi, d = 0)
}

barma.sim <- function(n, phi, seed, y.start = NULL){
  BARFIMA.sim(
    n = n,
    coefs = coeficientes(phi),
    y.start = y.start,
    error.scale = 1,
    complete = F,
    seed = seed
  )
}

barma.phi_estimado <- function(yt, alpha = 0, nu = 20, phi = 0.1){
  BARFIMA.fit(
    yt = yt,
    start = list(alpha = alpha, nu = nu, phi = phi),
    p = 1, # Odem do polinômio AR
    d = FALSE,
    error.scale = 1,
    report = F
  )$coefficients["phi"][[1]]
}

barma.residuos <- function(yt, phi_estimado){
  BARFIMA.extract(
    yt = yt,
    coefs = coeficientes(phi_estimado),
    llk = F,
    info = F,
    error.scale = 1
  )$error
}

# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func, debug = F, max = 5, ...) {
  FINALMENTE <- F
  contador <- 0
  
  while (!FINALMENTE) {
    contador <- contador + 1
    if (debug) print(paste("Tentativa:", contador))
    
    resultado <- tryCatch({
      func(...)
    }, error = function(err) {
      if (contador >= max) {
        stop("Deu ruim, número máximo de tentativas atingido")
      }
      
      if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
        # Ignora erro de extração de resíduos
        return(NULL)
      }
      
      stop(err)
    })
    
    if (!is.null(resultado)) {
      FINALMENTE <- T
    }
  }
  
  return(resultado)
}

Monte-Carlo

# Função auxiliar para gerar os dados
# Exemplo 1:
# ```r
# teste1 <- gerador_monte_carlo(
#   parametros = list(
#     lambda = seq(0.1, 0.4, by = 0.1)
#   ),
#   numero_de_execucoes = 2
# )
# 
# teste1
# ```
# 
# Exemplo 2:
# ```r
# teste2 <- gerador_monte_carlo(
#   parametros = expand.grid(
#     desvio_detectavel = seq(0.6, 1.8, by = 0.2),
#     intervalo_de_decisao = seq(4, 6, by = 1)
#   )
#   numero_de_execucoes = 2
# )
# 
# teste2
# ```
# 
# Total de execuções é:
# `execucoes = numero_de_execucoes * length(novos_phis) * nrow(parametros)`
# 
gerador_monte_carlo <- function(parametros = NULL, numero_de_execucoes = 200){
  # Verifica se `parametros` é um data frame, caso contrário converte para um data frame
  if (!is.null(parametros) && !is.data.frame(parametros)) {
    parametros <- as.data.frame(parametros)
  }
  
  novos_phis <- seq(0.2, 0.6, by = 0.1)
  
  # Expande a grade de parâmetros para cobrir todas as combinações
  result <- expand.grid(
    k = 1:numero_de_execucoes,
    h1_phi = novos_phis
  )
  
  if (!is.null(parametros)) {
    result <- merge(result, parametros, all = TRUE)
  }
  
  result %>%
    mutate(id = row_number()) %>%
    rowwise() %>%
    mutate(
      # H0
      ## Gera amostra de controle
      h0_amostras = list(barma.sim(
        n = nH0,
        phi = phi_parametro,
        seed = id
      )),
      ## Estima o valor de phi da amostra de controle
      h0_phi = barma.phi_estimado(h0_amostras),
      ## Calcula os resíduos da amostra de controle
      h0_residuos = list(barma.residuos(h0_amostras, h0_phi)),
      # H1
      ## Gera a amostra subsequente
      h1_controle = list(barma.sim(
        n = nH1,
        phi = h1_phi,
        # última observação da amostra de controle
        y.start = h0_amostras[nH0],
        seed = id + 1337E4
      )),
      ## Calcula os resíduos da amostra subsequente
      h1_residuos = list(so_quero_que_funcione(
        \() barma.residuos(h1_controle, h0_phi)
      ))
    )
}

Validação

Vamos verificar nossa implementação do gerador de Monte-Carlo.

teste_montecarlo <- cache_dados(
  "cache-teste-montecarlo",
  function() {
    gerador_monte_carlo(
      numero_de_execucoes = 200
    ) %>%
      select(h0_phi, h1_phi) %>%
      mutate(
        # Calcula a diferença entre os valores de phi
        # `phi_parametro`: valor de phi definido para a amostra de controle
        # `h0_phi`: valor de phi estimado para a amostra de controle. Espera-se que seja igual a `phi_parametro`
        diferenca = h0_phi - phi_parametro
      )
  }
)

ggplotly(
  teste_montecarlo %>%
    ggplot(aes(x = factor(1), y = diferenca)) +
    geom_boxplot() +
    geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
    annotate(
      geom = "text",
      x = 0.55,
      y = 0 + 0.01,
      label = "0.00",
      color = "red"
    ) +
    labs(
      x = "Valores de Φ₀",
      y = "Diferença entre Φ e Φ₀",
      title = paste0("Diferença entre Φ e Φ₀, para ", nrow(teste_montecarlo), " execuções")
    ) +
    theme.base
) %>%
  plotly.base

EWMA

O EWMA (Exponential Weighted Moving Average) é um método de controle de processos que utiliza uma média móvel ponderada exponencialmente para detectar mudanças no processo. Sendo definido por:

\[ z_i = \lambda y_i + (1 - \lambda) z_{i-1} \]

onde, \(z_i\) é a estatística de controle no instante \(i\), \(y_i\) é a observação no instante \(i\), \(\lambda\) é o fator de suavização e \(z_{i-1}\) é a estatística de controle no instante anterior.

O fator \(\lambda\) é uma constante definida no intervalo \((0, 1]\) e seu valor inicial (para \(i = 1\)) é definido como a média do processo, tal que \(z_0 = \mu_0\).

A estatística de controle, \(z_i\), é comparada com os limites de controle, \(\text{LCS}\) e \(\text{LCI}\), sendo é definido como:

\[ \text{LC} = \mu_0 \pm L \sigma \sqrt{\frac{\lambda}{2 - \lambda} \left[1 - \left(1 - \lambda\right)^{2i}\right]} \]

Aqui, utilizamos o pacote qcc para implementar o EWMA.

ewma_qcc <- function(amostra_inicial, lambda, amostra_subsequente) {
  ew <- ewma(
    amostra_inicial,
    lambda = lambda,
    nsigmas = quantil_teste,
    plot = F,
    newdata = amostra_subsequente
  )
  
  registros <- (nH0 + 1):(nH0 + nH1)
  
  ewma <- as.numeric(ew$y[registros])
  LI <- ew$limits[, 1][registros]
  LS <- ew$limits[, 2][registros]
  fora_de_controle <- ewma < LI | ewma > LS
  total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
  fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
  list(
    ewma = ewma,
    LI = LI,
    LS = LS,
    fora_de_controle = fora_de_controle,
    total_fora_de_controle = total_fora_de_controle,
    fracao_fora_de_controle = fracao_fora_de_controle
  )
}

Monte-Carlo

ewma_monte_carlo <- cache_dados(
  "cache-simulacao-ewma",
  function() {
    gerador_monte_carlo(
      parametros = list(
        lambda = seq(0.1, 0.4, by = 0.1)
      )
    ) %>%
      mutate(
        qcc = list(ewma_qcc(h0_residuos, lambda, h1_residuos)),
        ewma = list(qcc$ewma),
        LI = list(qcc$LI),
        LS = list(qcc$LS),
        fora_de_controle = list(qcc$fora_de_controle),
        total_fora_de_controle = qcc$total_fora_de_controle,
        fracao_fora_de_controle = qcc$fracao_fora_de_controle
      ) %>%
      select(-qcc)
  }
)

Cartas para λ = 0.2

Vamos analisar o EWMA para \(\lambda = 0.2\). Aqui, vamos considerar apenas a primeira execução.

ggplotly(
  ewma_monte_carlo %>%
    filter(k == 1 & lambda == 0.2) %>% # Apenas a primeira execução
    select(h1_phi, ewma, LI, LS, fora_de_controle, lambda) %>%
    rename(`Φ₁` = h1_phi) %>%
    mutate(
      observacao = list(1:nH1)
    ) %>%
    unnest(
      cols = c(`Φ₁`, observacao, ewma, LI, LS, fora_de_controle)
    ) %>%
    ggplot() +
    geom_line(aes(x = observacao, y = LI, color = "Limite Inferior"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = LS, color = "Limite Superior"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = ewma, color = "EWMA")) +
    geom_point(data = . %>% filter(fora_de_controle),
               aes(x = observacao, y = ewma, color = "Fora de controle"), size = 1.5, shape = 4) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
        "Limite Superior" = adjustcolor("red", alpha.f = 0.5),
        "EWMA" = adjustcolor("black", alpha.f = 0.8),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.4)
      )
    ) +
    labs(
      title = "EWMA para resíduos βAR(1), com λ = 0.2",
    ) +
    theme.base +
    theme(
      legend.position = "bottom",
      legend.title = element_blank()
    ) +
    facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
  plotly.base

Resumo

Resumo da fração de pontos fora de controle para diferentes valores de \(\lambda\) e \(\Phi_1\).

ewma_monte_carlo_resumo <- ewma_monte_carlo %>%
  group_by(lambda, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

datatable(
  ewma_monte_carlo_resumo %>%
   mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e λ",
  colnames = c("Valores de λ", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)

Gráficos

Análise gráfica das médias de FPFCs (Fração de Pontos Fora de Controle).

ggplotly(
  ewma_monte_carlo_resumo %>%
      ggplot(aes(x = h1_phi, y = mean, color = factor(lambda))) +
      geom_line() +
      geom_point() +
      labs(
        x = "Valores de Φ",
        y = "Fração de pontos fora de controle",
        color = "Valores de λ",
        title = "FPFCs por Φ e λ"
      ) +
      theme.base
) %>%
  plotly.base

E, a distribuiçõa dos PFCs.

ggplotly(
  ewma_monte_carlo %>%
    ggplot(aes(x = factor(h1_phi), y = fracao_fora_de_controle, fill = factor(lambda))) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      fill = "Valores de λ",
        title = "FPFCs por Φ e λ"
    ) +
    theme.base
  ) %>%
  plotly.base %>%
  layout(boxmode = 'group')

CUSUM

O CUSUM (Cumulative Sum) é um método de controle de processos que utiliza a soma acumulada dos resíduos para detectar mudanças no processo.

Neste método são definidas duas constantes, \(H\) e \(K\), que representam o tamanho da mudança que se deseja detectar e a sensibilidade do método, respectivamente. As estatísticas de controle, sendo duas, são definidas como:

\[ \begin{matrix} \text{C}^+_i & = & \max\left[0, x_i - (\mu_0 + K) + C^+_{i-1}\right] \\ \text{C}^-_i & = & \max\left[0, (\mu_0 - K) - x_i + C^-_{i-1}\right] \\ \end{matrix} \]

onde \(x_i\) é a observação no instante \(i\), \(\mu_0\) é a média do processo, \(C^+_0 = C^-_0 = 0\) e \(i = 1, 2, \ldots, n\).

A constante \(K\), chamada de valor de referência, geralmente, segundo Montgomery (2009), é definida como a metade entre o \(\mu_0\) alvo e o valor fora de controle de média \(\mu_1\) que queremos detectar.

Já a constante \(H\), chamada de limite de decisão, é definida como o valor que a estatística de controle deve atingir para que possamos detectar a mudança no processo. Para Montgomery (2009), um valor razoável para \(H\) é \(5\sigma\). Assim, se \(\text{C}^+_i \geq H\) ou \(\text{C}^-_i \leq -H\), então o processo está fora de controle.

cusum_qcc <- function(amostra_inicial_, desvio_detectavel = 1, intervalo_de_decisao = 5, amostra_subsequente = NULL) {
  arguments <- list(
    data = amostra_inicial_,
    se.shift = desvio_detectavel,
    decision.interval = intervalo_de_decisao,
    plot = F
  )
  
  if (!is.null(amostra_subsequente)) {
    arguments$newdata <- amostra_subsequente
  }
  
  cu <- do.call(cusum, arguments)
  
  registros <- if (is.null(amostra_subsequente)) {
    seq(1, nH0)
  } else {
    seq(nH0 + 1, nH0 + nH1)
  }
  pos <- cu[["pos"]][registros]
  neg <- cu[["neg"]][registros]
  fora_de_controle_pos <- pos > intervalo_de_decisao
  fora_de_controle_neg <- neg < -intervalo_de_decisao
  total_fora_de_controle <- sum(fora_de_controle_pos) + sum(fora_de_controle_neg)
  fracao_fora_de_controle <- total_fora_de_controle / length(registros)
  list(
    pos = pos,
    neg = neg,
    fora_de_controle_pos = fora_de_controle_pos,
    fora_de_controle_neg = fora_de_controle_neg,
    total_fora_de_controle = total_fora_de_controle,
    fracao_fora_de_controle = fracao_fora_de_controle
  )
}

Estimando o valor ótimo para K

Queremos que, sob \(H_1\), nas 100 amostras iniciais a probabilidade de um ponto fora de controle seja de, nominalmente, 5%. Portanto, vamos buscar um valor para K onde isso acontece.

desvio_detectavel_otimo_optimize <- function(amostra_inicial) {
  for (x in seq(0, 3, by = 0.1)) {
    fora <- cusum_qcc(amostra_inicial, desvio_detectavel = x)$total_fora_de_controle
    if (fora == 5) {
      return(x)
    }
  }
  return(NA)
}

cusum_k_otimo <- cache_dados(
  "cache-cusum-k-otimo",
  \() {
    data.frame(
      id = 1:1000
    ) %>%
    rowwise() %>%
    mutate(
      h0_amostra = list(barma.sim(
        n = nH0,
        phi = phi_parametro,
        seed = id
      )),
      h0_phi = list(barma.phi_estimado(h0_amostra)),
      h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
      minimo = desvio_detectavel_otimo_optimize(h0_residuos)
    )
  }
)

ggplotly(
  cusum_k_otimo %>%
    filter(!is.na(minimo)) %>%
    ggplot(
      aes(x = 0, y = minimo, group = 0, fill = "red")
    ) +
    geom_boxplot() +
    labs(
      x = element_blank(),
      y = "Valor ótimo de K",
      title = "Valor ótimo de K por simulação"
    ) +
    theme.base +
    theme.no_legend
) %>%
  plotly.base

valor_otimo_k <- mean(cusum_k_otimo$minimo, na.rm = TRUE)

kableExtra::kbl(
  cusum_k_otimo %>%
    filter(!is.na(minimo)) %>%
    group_by() %>%
    summarise(
      "Média" = mean(minimo),
      "Desvio padrão" = sd(minimo),
      "Amostras" = n(),
    ),
  caption = "Valor ótimo para K",
  align = 'c',
  digits = 3
)

Valor ótimo para K

Média	Desvio padrão	Amostras
0.507	0.193	260

Monte-Carlo

cumsum_k_monte_carlo <- cache_dados(
  "cache-simulacao-cusum-k",
  function() {
    gerador_monte_carlo(
      parametros = list(
        desvio_detectavel = seq(0.1, 1.0, by = 0.1)
      )
    ) %>%
      mutate(
        qcc = list(cusum_qcc(h0_residuos, desvio_detectavel = desvio_detectavel, amostra_subsequente = h1_residuos)),
        pos = list(qcc$pos),
        neg = list(qcc$neg),
        fora_de_controle_pos = list(qcc$fora_de_controle_pos),
        fora_de_controle_neg = list(qcc$fora_de_controle_neg),
        total_fora_de_controle = qcc$total_fora_de_controle,
        fracao_fora_de_controle = qcc$fracao_fora_de_controle
      ) %>%
      select(-qcc)
  }
)

Cartas para K ótimo

Vamos analisar o CUSUM para o valor ótimo de K encontrado.

cusum_para_k_otimo <- cumsum_k_monte_carlo %>%
  filter(k == 1 & desvio_detectavel == round(valor_otimo_k, 1)) %>%
  rename(`Φ₁` = h1_phi) %>%
  mutate(
    observacao = list(1:nH1)
  )

ggplotly(
  cusum_para_k_otimo %>%
    unnest(
      cols = c(`Φ₁`, observacao, pos, neg, h1_residuos, fora_de_controle_pos, fora_de_controle_neg)
    ) %>%
    ggplot() +
    geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
    geom_hline(aes(yintercept = -5, color = "Limite de decisão"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = pos, color = "N+")) +
    geom_line(aes(x = observacao, y = neg, color = "N-")) +
    geom_point(data = . %>% filter(fora_de_controle_pos),
               aes(x = observacao, y = pos, color = "Fora de controle"), size = 1.5, shape = 4) +
    geom_point(data = . %>% filter(fora_de_controle_neg),
               aes(x = observacao, y = neg, color = "Fora de controle"), size = 1.5, shape = 4) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "N+" = adjustcolor("blue", alpha.f = 0.6),
        "N-" = adjustcolor("darkgreen", alpha.f = 0.6),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.3),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.4),
        "Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
      )
    ) +
    labs(
      title = "CUSUM para resíduos βAR(1)",
    ) +
    theme.base +
    theme(
      legend.position = "bottom",
      legend.title = element_blank()
    ) +
    facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
  plotly.base

Resumo

Em média, diminuir o K aumenta a sensibilidade do CUSUM para detectar mudanças no processo.

cumsum_k_monte_carlo_resumo <- cumsum_k_monte_carlo %>%
  group_by(desvio_detectavel, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

datatable(
  cumsum_k_monte_carlo_resumo %>%
   mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e K, com H = 5",
  colnames = c("Valores de K", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)

Gráficos

Análise gráfica das médias de FPFCs.

ggplotly(
  cumsum_k_monte_carlo_resumo %>%
      ggplot(aes(x = h1_phi, y = mean, color = factor(desvio_detectavel))) +
      geom_line() +
      geom_point() +
      labs(
        x = "Valores de Φ₁",
        y = "Fração de pontos fora de controle",
        color = "Valores de K",
        title = "FPFCs por Φ₁ e K com H = 5"
      ) +
      theme.base
) %>%
  plotly.base

E, a distribuição dos PFCs.

ggplotly(
  cumsum_k_monte_carlo %>%
    ggplot(aes(x = factor(h1_phi), y = fracao_fora_de_controle, fill = factor(desvio_detectavel))) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      fill = "Valores de K",
        title = "FPFCs por Φ₁ e K com H = 5"
    ) +
    theme.base
  ) %>%
  plotly.base %>%
  layout(boxmode = 'group')

Estimando o valor ótimo para H

Semelhante ao valor ótimo para K, queremos que, sob \(H_1\), nas 100 amostras iniciais a probabilidade de um ponto fora de controle seja de, nominalmente, 5%. Portanto, vamos buscar um valor para H onde isso acontece.

intervalo_de_decisao_otimo_optimize <- function(amostra_inicial) {
  for (x in seq(1, 6, by = 1)) {
    fora <- cusum_qcc(amostra_inicial, intervalo_de_decisao = x)$total_fora_de_controle
    if (fora == 5) {
      return(x)
    }
  }
  return(NA)
}

cusum_h_otimo <- cache_dados(
  "cache-cusum-h-otimo",
  \() {
    data.frame(
      id = 1:1000
    ) %>%
    rowwise() %>%
    mutate(
      h0_amostra = list(barma.sim(
        n = nH0,
        phi = phi_parametro,
        seed = id
      )),
      h0_phi = list(barma.phi_estimado(h0_amostra)),
      h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
      minimo = intervalo_de_decisao_otimo_optimize(h0_residuos)
    )
  }
)

ggplotly(
  cusum_h_otimo %>%
    filter(!is.na(minimo)) %>%
    ggplot(
      aes(x = 0, y = minimo, group = 0, fill = "red")
    ) +
    geom_boxplot() +
    labs(
      x = element_blank(),
      y = "Valor ótimo de K",
      title = "Valor ótimo de K por simulação"
    ) +
    theme.base +
    theme.no_legend
) %>%
  plotly.base

valor_otimo_h <- mean(cusum_h_otimo$minimo, na.rm = TRUE)

kableExtra::kbl(
  cusum_h_otimo %>%
    filter(!is.na(minimo)) %>%
    group_by() %>%
    summarise(
      "Média" = mean(minimo),
      "Desvio padrão" = sd(minimo),
      "Amostras" = n(),
    ),
  caption = "Valor ótimo para H",
  align = 'c',
  digits = 3
)

Valor ótimo para H

Média	Desvio padrão	Amostras
3.16	0.712	125

Monte-Carlo

cumsum_h_monte_carlo <- cache_dados(
  "cache-simulacao-cusum-h",
  function() {
    gerador_monte_carlo(
      parametros = list(
        intervalo_de_decisao = seq(1, 6, by = 1)
      )
    ) %>%
      mutate(
        qcc = list(cusum_qcc(h0_residuos, intervalo_de_decisao = intervalo_de_decisao, amostra_subsequente = h1_residuos)),
        pos = list(qcc$pos),
        neg = list(qcc$neg),
        fora_de_controle_pos = list(qcc$fora_de_controle_pos),
        fora_de_controle_neg = list(qcc$fora_de_controle_neg),
        total_fora_de_controle = qcc$total_fora_de_controle,
        fracao_fora_de_controle = qcc$fracao_fora_de_controle
      ) %>%
      select(-qcc)
  }
)

Cartas para K ótimo

Vamos analisar o CUSUM para o valor ótimo de K encontrado.

cusum_para_h_otimo <- cumsum_h_monte_carlo %>%
  filter(k == 1 & intervalo_de_decisao == round(valor_otimo_h)) %>%
  rename(`Φ₁` = h1_phi) %>%
  mutate(
    observacao = list(1:nH1)
  )

ggplotly(
  cusum_para_h_otimo %>%
    unnest(
      cols = c(`Φ₁`, observacao, pos, neg, h1_residuos, fora_de_controle_pos, fora_de_controle_neg)
    ) %>%
    ggplot() +
    geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
    geom_hline(aes(yintercept = -5, color = "Limite de decisão"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = pos, color = "N+")) +
    geom_line(aes(x = observacao, y = neg, color = "N-")) +
    geom_point(data = . %>% filter(fora_de_controle_pos),
               aes(x = observacao, y = pos, color = "Fora de controle"), size = 1.5, shape = 4) +
    geom_point(data = . %>% filter(fora_de_controle_neg),
               aes(x = observacao, y = neg, color = "Fora de controle"), size = 1.5, shape = 4) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "N+" = adjustcolor("blue", alpha.f = 0.6),
        "N-" = adjustcolor("darkgreen", alpha.f = 0.6),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.3),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.4),
        "Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
      )
    ) +
    labs(
      title = "CUSUM para resíduos βAR(1)",
    ) +
    theme.base +
    theme(
      legend.position = "bottom",
      legend.title = element_blank()
    ) +
    facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
  plotly.base

Resumo

Em média, diminuir o H aumenta a sensibilidade do CUSUM para detectar mudanças no processo.

cumsum_h_monte_carlo_resumo <- cumsum_h_monte_carlo %>%
  group_by(intervalo_de_decisao, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

datatable(
  cumsum_h_monte_carlo_resumo %>%
   mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e H, com K = 1",
  colnames = c("Valores de H", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)

Gráficos

Análise gráfica das médias de FPFCs.

ggplotly(
  cumsum_h_monte_carlo_resumo %>%
      ggplot(aes(x = h1_phi, y = mean, color = factor(intervalo_de_decisao))) +
      geom_line() +
      geom_point() +
      labs(
        x = "Valores de Φ₁",
        y = "Fração de pontos fora de controle",
        color = "Valores de H",
        title = "FPFCs por Φ₁ e H, com K = 1",
      ) +
      theme.base
) %>%
  plotly.base

E, a distribuição dos PFCs.

ggplotly(
  cumsum_h_monte_carlo %>%
    ggplot(aes(x = factor(h1_phi), y = fracao_fora_de_controle, fill = factor(intervalo_de_decisao))) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      fill = "Valores de K",
        title = "FPFCs por Φ₁ e K com H = 5"
    ) +
    theme.base
  ) %>%
  plotly.base %>%
  layout(boxmode = 'group')

EWMA x CUSUM(K/H)

Vamos comparar os algoritmos EWMA e CUSUM(K) e CUSUM(H) para diferentes valores de \(\Phi_1\).

Aqui chamamos de CUSUM(K) as estatísticas onde o \(H\) for mantido constante e o \(K\) é variado, e o inverso para CUSUM(H).

estatisticas_resumo_combinadas <- bind_rows(
  cumsum_k_monte_carlo_resumo %>%
    filter(desvio_detectavel == 0.9) %>%
    rename(parametro = desvio_detectavel) %>% mutate(algoritmo = "CUSUM(K)"),
  cumsum_h_monte_carlo_resumo %>%
    filter(intervalo_de_decisao == 4) %>%
    rename(parametro = intervalo_de_decisao) %>% mutate(algoritmo = "CUSUM(H)"),
  ewma_monte_carlo_resumo %>%
    filter(lambda == 0.4) %>%
    rename(parametro = lambda) %>% mutate(algoritmo = "EWMA")
) %>%
  select(algoritmo, parametro, h1_phi, mean) %>%
  mutate(
    parametro = as.factor(parametro)
  )

ggplotly(
  estatisticas_resumo_combinadas %>%
    ggplot() +
    geom_line(
      aes(x = h1_phi, y = mean, color = parametro, linetype = algoritmo)
    ) +
    geom_point(
      aes(x = h1_phi, y = mean, color = parametro, shape = algoritmo)
    ) +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros",
      shape = "Algoritmo"
    ) +
    theme.base
) %>%
  plotly.base

estatisticas_combinadas <- bind_rows(
  cumsum_k_monte_carlo %>%
    filter(desvio_detectavel == 0.9) %>%
    rename(parametro = desvio_detectavel) %>% mutate(algoritmo = "CUSUM(K)"),
  cumsum_h_monte_carlo %>%
    filter(intervalo_de_decisao == 4) %>%
    rename(parametro = intervalo_de_decisao) %>% mutate(algoritmo = "CUSUM(H)"),
  ewma_monte_carlo %>%
    filter(lambda == 0.4) %>%
    rename(parametro = lambda) %>% mutate(algoritmo = "EWMA")
) %>%
  select(algoritmo, parametro, h1_phi, fracao_fora_de_controle) %>%
  mutate(
    parametro = as.factor(parametro),
    h1_phi = as.factor(h1_phi)
  )

ggplotly(
  estatisticas_combinadas %>%
    ggplot() +
    geom_boxplot(
      aes(x = h1_phi, y = fracao_fora_de_controle, fill = parametro, shape = algoritmo),
    ) +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      fill = "Valores dos parâmetros",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros"
    ) +
    facet_wrap(~algoritmo) +
    theme.base
) %>%
  plotly.base %>%
  layout(boxmode = 'group')

Análise

A partir da análise da média dos pontos fora de controle, observa-se que o CUSUM(K), em média, possui poder semelhante ao CUSUM(H), enquanto que o EWMA é inferior aos dois.

Já a análise gráfica dos boxplots mostra que o CUSUM(K) possui uma variabilidade, em geral, maior que o CUSUM(H), enquanto que o EWMA possui, em geral, a menor variabilidade entre os algoritmos.

E se… juntarmos CUSUM(K) e CUSUM(H)?

Vamos verificar se existe alguma combinação ótima entre os valores de K e H para detectar pontos fora de controle.

comb_cusum <- cache_dados(
  "cache-simulacao-cusum-k-h",
  function() {
    gerador_monte_carlo(
      parametros = expand.grid(
        desvio_detectavel = seq(0.6, 1.8, by = 0.2),
        intervalo_de_decisao = seq(3, 6, by = 1)
      )
    ) %>%
      mutate(
        qcc = list(cusum_qcc(h0_residuos,
                             desvio_detectavel = desvio_detectavel,
                             intervalo_de_decisao = intervalo_de_decisao,
                             amostra_subsequente = h1_residuos)),
        pos = list(qcc$pos),
        neg = list(qcc$neg),
        fora_de_controle_pos = list(qcc$fora_de_controle_pos),
        fora_de_controle_neg = list(qcc$fora_de_controle_neg),
        total_fora_de_controle = qcc$total_fora_de_controle,
        fracao_fora_de_controle = qcc$fracao_fora_de_controle
      ) %>%
      select(-qcc)
  }
)

Resultados

comb_cusum_resumo <- comb_cusum %>%
  group_by(desvio_detectavel, intervalo_de_decisao, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  ) %>%
  mutate(
    desvio_detectavel = as.factor(desvio_detectavel),
    intervalo_de_decisao = as.factor(intervalo_de_decisao),
  )

datatable(
  comb_cusum_resumo %>%
   mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle",
  colnames = c("Valores de K", "Valores de H", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)

ggplotly(
  comb_cusum_resumo %>%
    group_by(desvio_detectavel, intervalo_de_decisao) %>%
    summarise(
      h1_phi = list(h1_phi),
      mean = list(mean),
      .groups = "drop"
    ) %>%
    filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
    unnest(cols = c(h1_phi, mean)) %>%
    ggplot(aes(x = h1_phi, y = mean, color = desvio_detectavel, linetype = intervalo_de_decisao)) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      color = "Valores de K",
      linetype = "Valores de H",
      title = "FPFCs por Φ₁, K e H"
    ) +
    theme.base
) %>%
  plotly.base

De uma forma geral, observa-se que aumentando K e H diminui a sensibilidade do CUSUM para detectar pontos fora de controle.

E, analisando juntamente com as outras execuções, temos

estatisticas_combinadas_resumo_kh <- bind_rows(
  estatisticas_resumo_combinadas,
  comb_cusum_resumo %>%
    filter(desvio_detectavel == 0.8 & intervalo_de_decisao == 5) %>%
    mutate(parametro = as.factor(paste0(desvio_detectavel, ";", intervalo_de_decisao))) %>%
    mutate(algoritmo = "CUSUM(K/H)") %>%
    select(algoritmo, parametro, h1_phi, mean)
)

ggplotly(
  estatisticas_combinadas_resumo_kh %>%
    ggplot() +
    geom_line(
      aes(x = h1_phi, y = mean, color = parametro, linetype = algoritmo)
    ) +
    geom_point(
      aes(x = h1_phi, y = mean, color = parametro, shape = algoritmo)
    ) +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros",
      shape = "Algoritmo"
    ) +
    theme.base
) %>%
  plotly.base

estatisticas_combinadas_kh <- bind_rows(
  estatisticas_combinadas,
  comb_cusum %>%
    filter(desvio_detectavel == 0.8 & intervalo_de_decisao == 5) %>%
    mutate(parametro = paste0(desvio_detectavel, ";", intervalo_de_decisao)) %>%
    mutate(algoritmo = "CUSUM(K/H)") %>%
    select(algoritmo, parametro, h1_phi, fracao_fora_de_controle) %>%
    mutate(
      parametro = as.factor(parametro),
      h1_phi = as.factor(h1_phi)
    )
)

ggplotly(
  estatisticas_combinadas_kh %>%
    ggplot() +
    geom_boxplot(
      aes(x = h1_phi, y = fracao_fora_de_controle, fill = parametro, shape = algoritmo),
    ) +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      fill = "Valores dos parâmetros",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros"
    ) +
    facet_wrap(~algoritmo) +
    theme.base
) %>%
  plotly.base %>%
  layout(boxmode = 'group')

De uma forma geral, observa-se um leve aumento no poder quando combinamos H = 5 e K = 0.8. Além diss, houve um leve aumento na variabilidade.

E se… combinarmos os dois algoritmos?

Vamos combinar os dois algoritmos e ver o que acontece. A ideia é que, se um dos algoritmos detectar um ponto fora de controle, então o processo está fora de controle.

comb_monte_carlo <- cache_dados(
  "cache-simulacao-ewma-cusum",
  function() {
    gerador_monte_carlo(
      parametros = expand.grid(
        desvio_detectavel = seq(1.6, 2.0, by = 0.2),
        intervalo_de_decisao = seq(4, 6, by = 1),
        lambda = seq(0.6, 1.0, by = 0.2)
      )
    ) %>%
      mutate(
        "(K;H;λ)" = factor(paste0(desvio_detectavel, ";", intervalo_de_decisao, ";", lambda)),
        ewma = list(ewma_qcc(h0_residuos, lambda, h1_residuos)),
        cumsum = list(cusum_qcc(h0_residuos,
                                desvio_detectavel = desvio_detectavel,
                                intervalo_de_decisao = intervalo_de_decisao,
                                amostra_subsequente = h1_residuos)),
        fora_de_controle = list(
          ewma$fora_de_controle | (cumsum$fora_de_controle_pos | cumsum$fora_de_controle_neg)
        ),
        fracao_fora_de_controle = sum(fora_de_controle) / nH1
      ) %>%
      select(-ewma, -cumsum)
  }
)

Resultados

comb_monte_carlo_resumo <- comb_monte_carlo %>%
  group_by(`(K;H;λ)`, h1_phi, desvio_detectavel, lambda) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

datatable(
  comb_monte_carlo_resumo %>%
    select(-`(K;H;λ)`) %>%
    mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e (K;H;λ)",
  colnames = c("Valores de Φ₁", "K", "H", "λ", "Média", "Mínimo", "Máximo")
)

ggplotly(
  comb_monte_carlo_resumo %>%
    # get `(K;H;λ)` where h1_phi == 0.2 and mean >= 0.05 & mean < 0.06
    group_by(`(K;H;λ)`) %>%
    summarise(
      h1_phi = list(h1_phi),
      mean = list(mean),
      .groups = "drop"
    ) %>%
    filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
    unnest(cols = c(h1_phi, mean)) %>%
    ggplot(aes(x = h1_phi, y = mean, color = `(K;H;λ)`)) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      color = "Constantes (K;H;λ)",
      title = "FPFC por Φ₁ e (K;H;λ)"
    ) +
    theme.base
) %>%
  plotly.base

A melhor combinação encontrada

estatisticas_combinadas_resumo_ew_sum <- bind_rows(
  estatisticas_combinadas_resumo_kh,
  comb_monte_carlo_resumo %>%
    filter(`(K;H;λ)` == "1.8;6;0.8") %>%
    mutate(parametro = as.factor(`(K;H;λ)`)) %>%
    mutate(algoritmo = "EW-SUM") %>%
    select(algoritmo, parametro, h1_phi, mean)
)

ggplotly(
  estatisticas_combinadas_resumo_ew_sum %>%
    ggplot() +
    geom_line(
      aes(x = h1_phi, y = mean, color = parametro, linetype = algoritmo)
    ) +
    geom_point(
      aes(x = h1_phi, y = mean, color = parametro, shape = algoritmo)
    ) +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros",
      shape = "Algoritmo"
    ) +
    theme.base
) %>%
  plotly.base

Análise

De uma forma geral, houve uma piora na detecção de pontos fora de controle ao combinar os dois algoritmos. Isso pode ser explicado pelo fato de que, a sensibilidade do EWMA é maior que a do CUSUM, o que pode ter levado a um aumento de falsos positivos quando o processo permanece sob controle.